

%% Load WRF output ***instantaneous variable***
%% Take U wind as an example

tic

% Loc
Loc = '.\';

% sec
Sce = { '1995', '2018' };

% years
time_year = 2018 : 2020;

% interp levels in WRF
interp_levels = [500,400,300,200];

% select load level
select_PL_level = [500,400,300,200];
flag_level = find( sum( interp_levels==select_PL_level, 2 ) );

% Load
W_single_PL_U_aft = zeros( [size(W_lon), length( select_PL_level ), 365, length(time_year) ] );	% UNEP
W_single_PL_U_dif = zeros( [size(W_lon), length( select_PL_level ), 365, length(time_year) ] );	% CTRL

for y = 1:size(time_year,2)
   
    for m = 1 : 12

        md = eomday(2018, m);

        for d = 1:md

            n = sum( eomday(2018, [1:m-1]) ) + d;

            % read
            file_bef = [ Loc, 'a_', num2str( time_year(y) ), '_', Sce{ 1 }, '\wrfpress_d02_', ...
                    num2str(time_year(y)), '-', num2str(m,'%02d'), '-', num2str(d,'%02d'), '_00_00_00.nc' ];
    
            file_aft = [ Loc, 'a_', num2str( time_year(y) ), '_', Sce{ 2 }, '\wrfpress_d02_', ...
                    num2str(time_year(y)), '-', num2str(m,'%02d'), '-', num2str(d,'%02d'), '_00_00_00.nc' ];
    
            PL_U_bef = rot90(  ncread(file_bef, 'Umet_PL')  );
            PL_U_bef = PL_U_bef(:,:,flag_level);
            PL_U_aft = rot90(  ncread(file_aft, 'Umet_PL')  );
            PL_U_aft = PL_U_aft(:,:,flag_level);
            
            W_single_PL_U_bef(:,:,:,n,y) = PL_U_bef;
            W_single_PL_U_aft(:,:,:,n,y) = PL_U_aft;

        end

        disp( [num2str(time_year(y)), '-', num2str(m,'%02d')] )
        
    end

end

clear PL_U_bef PL_U_aft 

toc 






